Structural and functional networks in complex systems with delay 
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Functional networks of complex systems are obtained from the analysis of the temporal activity 
of their components, and are often used to infer their unknown underlying connectivity. We obtain 
the equations relating topology and function in a system of diffusively delay-coupled elements in 
complex networks. We solve exactly the resulting equations in motifs (directed structures of three 
nodes), and in directed networks. The mean- field solution for directed uncorrelated networks shows 
that the clusterization of the activity is dominated by the in-degree of the nodes, and that the 
locking frequency decreases with increasing average degree. We find that the exponent of a power 
law degree distribution of the structural topology, 7, is related to the exponent of the associated 
functional network as a = (2 — 7)"^, for 7 < 2. 

PACS numbers: 89.75.Fb, 05.40.-a, 05.65.+b, 89. 75. He 
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I. INTRODUCTION 

Collective phenomena in populations of interacting el- 
ements is a subject of intense study in physical, biolog- 
ical, chemical, and social systems [l|-y]- In many cases, 
the emergence of patches of coherent behavior is the main 
observable we have of the underlying dynamics and inter- 
action of their constituents. This is the case, for example, 
in gene expression, measured as DNA levels in microar- 
rays or electrophysiological activity in the brain @, 

measured through multi-unit extracellular electrode. In 
many occasions, the coordination is not global but local, 
and the observation reveals clusters of elements dynam- 
ically correlated or, generally speaking, synchronized 0- 
[To| . The resulting networks of coordinated activity are 
usually called functional networks of the system [ij, [l^] ■ 
The analysis of these networks, from the physicists per- 
spective, commonly focuses on network synchrony in the 
absence of time delays. However, delays are common 
in neural networks [13[, and many other biological 
[Tg} and social systems where the interaction between ele- 
ments involves the propagation through a communication 
channel. The consideration of these delays is of utmost 
importance [13] ■ Recently, it has been analytically, and 
experimentally shown that zero time-lag synchronization 
is feasible over two distant (delayed) interacting oscilla- 
tors when a third oscillator is placed in between of them 
[isl [l9| . A recent study sheds light along these lines by 
studying the synchronization of networks of chaotic units 
with time-delayed couplings using the formalism of the 
master stability function j20j] . 

There has been a big effort from the scientific commu- 
nity to infer the complex network of interactions between 
elements from the functional network [2l|-|28ll. Here we 



investigate on this inference from a fundamental physical 
perspective. We analyze the functional network resulting 
from the simplest dynamical system with delay present- 
ing a synchronous dynamics, on a given topology, and 
relate topology and functionality. Given the simplicity 
of the model, we obtain the exact solution, develop a 
statistical mean-field theory approximation and find the 
relation between the degree distribution of the topologi- 
cal network and the associated functional network. 



II. THE MODEL 

First, we develop the analytical aspects of the prob- 
lem. Let's start considering a set of N elements coupled 
diffusively with delay 



(f>i{t) = UJi + 



(1) 



where Uij are the components of the adjacency matrix A, 
that is, ttij = 1 if element j influences i with a delay Tij, 
and e is the coupling strength. Without loss of general- 
ity we can rescale the coupling strength to e = 1. Equa- 
tion ll]) represents a system of N elements in a network, 
which move at constant speed and that adjust their local 
position to match that of their neighbors; the communi- 
cation between pair of nodes is not instantaneous but it is 
characterized by some delay. An alternative interpreta- 
tion of Eq. ([T]) corresponds to the linearization of many 
non-linear interaction models, including the Kuramoto 
model (2^ , as long as the phase differences remain small 
enough. As we will demonstrate later, if the network can 
be reached from at least one node, Eq. ([1]) presents a 
unique phase locked solution of the form 



(2) 
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where fl is the locking frequency and 9i the initial phase 
of element i. Substituting in Eq. ([T]), we obtain a set of 
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N linear equations that can be written in matrix form as 



a; - r2 (1 + T) = 



(3) 



where L is the Laplacian matrix defined as Lij — 
ail, ki.in = J2j '^ij tti^ in-degree of node i, 



Sij is the Kronecker delta, 1 is a vector of 1 's and T is a 
vector of components Ti = atjTij the total delay af- 
fecting each node, u: and 6 are the frequency vector and 
the phase vector with components uJi and 9i respectively. 

In general, the Laplacian matrix L is asymmetric, how- 
ever, as the sum of its rows is zero, it admits a left- 
eigenvector c = (ci, C2, Cat) with eigenvalue 0, that is, 
cL = 0. The left eigenvector c is unique as long as all 
the nodes of the network can be reached by at least one 
node [30[. Left-muhiplying Eq. dS]) by c we obtain the 
locking frequency 



1 + (T) ' 



(4) 



where (x) = CiXt and c is normalized, J2i = 1- The 
phases are now given by 



u:, - (u;) + (c^,(T) - (u;)T,) 
1 + (T) 



(L0) 



(5) 



For undirected networks, the left-eigenvector c = 
{1/N, 1/N, is unique if it contains a single com- 
ponent and the brackets in Eq. arc unweighted aver- 
ages. 

Thus, if all the nodes of the network can be reached by 
at least one node, then Eq. ([T]) has a unique phase locked 
solution given by Eqs. (|4]) and ([5|). These equations point 
out the relationship between the topology of the network, 
the distribution of delays, the locking frequency, and the 
state of the elements. 



A. Perfect synchronization 

The topology, delays and frequency can be combined 
to achieve the same state for each node. The condition 
for fully clustered solutions {9i = 9j, implies 



(6) 



From Eq. ([6]) we see that in absence of delays {Ti — 0) all 
the elements must have the same natural frequency, and 
reversely, if all nodes have the same natural frequency 
(wi = uj), the total delay affecting each node must be 
the same. It is straightforward to prove that Eq. ^ 
is satisfied for degree regular networks of identical ele- 
ments with equal delays, i.e., oJi = oj and Ti = r. In 
this case, the frequency of the phase synchronized state 
is = TTrfc' general, it is always possible to choose 
the frequencies, the topology and the delays such that 
perfect coordinated activity is reached. In this case, the 
functional network, that is, the network formed connect- 
ing those nodes displaying correlated activity, will be a 
fully connected network despite the sparse connectivity 
of the underlying interaction network. 



B. Motifs 

Beyond the above perfectly synchronized solutions, we 
will extend our analysis to directed networks paying at- 
tention to the clusterization of the activity with different 
values of the phases. The simplest possible case corre- 
sponds to graphs of three nodes, motifs HH, HI] . The 
interest in motifs comes from previous studies showing 
the impact of motifs synchronization in absence of de- 
lays as building blocks of larger synchronized structures 
|33| . Using Eq. we find the locked solution for each of 
the twelve different motifs of three elements with directed 
couplings. For the sake of clarity, we assume each element 
has the same frequency LOi = ui, and delay Tij = t, Vi, j. 
Solving Eq. ([3]) we obtain for every motif configuration 
the normalized oscillation frequency r~^ = fl/ui and the 



phase differences 



,) ~ fir Ay. The twelve different 



motifs are classified in five different functional networks 
(see Fig. [T|) . This result points out the impossibility of 
deriving the motif topology solely from the information 
of the functional networks due to the degeneracy shown 

a. 



III. HETEROGENEOUS MEAN FIELD 
APPROACH 

Beyond the formal exact solution presented in Eqs. 
we want to gain insight on the class of uncor- 
related directed networks. First we start considering 
heterogeneous directed networks, specified by their de- 
gree distribution /'(k), where k = (fcin, fcout), and by 
the conditional probability P(k'|k) that a node of de- 
gree k is connected to a node of degree k'. Normal- 
ization conditions J^w^O^) ~ ^ ^^'^ J2w ~ ^ 
must be fulfilled. The degree detailed balance condition 
A;outP(k)Pout(k'|k) = A:i'„P(k')Pin(k|k') (where Pi„(k'|k) 
[Pout(k'|k)] measures the probability to reach a vertex 
of degree k' leaving from a vertex of degree k using an 
incoming [outgoing] edge of the source vertex) ensures 
that the network is closed and that (kin) = (fcout)- We 
resort on the heterogeneous mean-field approach, coarse- 
graining the dynamics to classes of nodes of the same 
degree k. Thus, we define the phase density of nodes 
of degree k as 



(7) 



where iVk = P(k)N is the expected number of nodes 
with degree k. Here we have made use of K to denote 
the set of nodes with degree k. Similarly we define the 
frequency density 



1 



i£K 



(8) 



This notation allows to group the sums by the degrees of 
the nodes. For instance, if the degree of node ? is k^ = k 
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FIG. 1: Functional network (green) for twelve structural motifs of three elements. In the functional solution the black solid 
line represents a zero phase difference between the oscillators while the gray dashed lines stand for a non-zero phase difference. 
For the structural, solid black lines stand for interaction links and directionality is indicated by an arrow. Parameters: Aij = 
(Oi - 9j)/flT, and r = oj/fi. 



then 



a,j(j)j = kin Pin(k'|k)$k' 



(9) 



k' 



For identical elements uji = to and Tij — t, the time 
evolution of the phase density of the class of nodes of 
degree k, ^k(t), can be rewritten from Eq. ([Ij as 

= W^k + fci„^Pi„(k'|k)($k'(i-T)-$kW) • (10) 

k' 

For uncorrelated networks Pin(k'|k) = and 
with the ansatz of locked solutions $k = + ©k, we 
obtain 

n = w^- hMr + E fc;ut^(k')(ek' - Ok) . (11) 

Summing over all degrees wc find 

n{k)T 1 



ew = 



n = 



I 



{feout/fc,„) 



(12) 
(13) 



being a an arbitrary constant. For undirected networks 
hn = fcout, thus we recover Eq. (jH) for the locking fre- 
quency where (fc)r = (T), and 



n{k}T 

Bfe = — ; h a 



This indicates that whether two nodes show a similar 
phase depends on their degree difference in an uncorre- 
lated network. It also shows that low-degree nodes are 
ahead of high-degree nodes. At least in this limit, the 
precise shape of the degree distribution is not playing an 
important role, as only the average degree (k) enters into 
the equation. Obviously, this dependence of the degree 
is reminiscent of our hypothesis of a mean-field coarse- 
grained by degree, however it is not trivial that this ap- 
proximation will hold for the actual dynamics (Eq. p])). 

Further more the distribution of the phases in cor- 
related networks as for example in the C. elegans neu- 
ral network also shows a good agreement with Eq. (|14p . 
The neuronal network connectivity of the C. elegans can 
be represented as a weighted adjacency matrix of 275 
nonpharyngeal neurons, out of a total of 302 neurons 



(14) 



( |http: / / www.wormatlas.org/ 1 . We assume that the ner- 
vous system of the C. elegans can be modeled as a net- 
work, where nodes represent the center of the cell bod- 
ies, and the links represent synapses. The heterogeneous 
mean field formalism describes the relationship between 
dynamics and topology in uncorrelated networks. Such 
relationship can be illustrated in a real (correlated) net- 
work analyzing the dynamics of Eq. (1) using the con- 
nectivity of the neural system of the C. elegans (Fig. [2]). 
When comparing the exact solution in the directed neu- 
ral network with the analytical solution we observe that 
it captures the dependence on the in-degree and gives an 
excellent solution for the rewired directed network. Thus 
the in-degree of a neuron gives a good first approximation 
to the real state of the neuron although the precise wiring 
details are very important to know its exact value. Sim- 
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FIG. 2: (Top row) Phase clusterization in the directed C. 
elegans neural network. We plot the phases versus the in- 
degree kin for neurons in the C. elegans (open circles) and 
rewired networks (filled symbols) (averaged over f 00 realiza- 
tions of the rewiring algorithm) keeping the same (/cin, /i^out) 
for each neuron. In the rewired networks the phase is well 
approximated by the relation 0{k) — b/k + a. The phases 
are obtained after integration of Eq. (1) with uji = 1 and 
Tij = 0.1. (Bottom row) Adjacency matrix of the neuronal 
connections (a), and of the functional network (b). In both 
cases the neurons are ordered according to the ranking of their 
phases obtained from the dynamical system given by Eq. (1) 
with uui = 1 and Tij =0.1. 
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FIG. 3: Functional degree distribution P{q) for a network 
of 10^ nodes with structural degree distribution P{k) ~ k~'^ 
and 7 = 1.5. Dashed line corresponds to q~^. Inset: De- 
pendence of (circles) and a (squares) on the structural 
degree distribution exponent 7. Solid black line represents the 
theoretical prediction. The functional networks are obtained 
after integration of Eq. ([1} with uji = 1 and Tij = 0.1, and 
A = 10-^ 



agreement with the analytical prediction for classes of 
nodes with degree k in uncorrelated networks. In Fig. [3] 
we compare the values of the exponent a and /3 obtained 
after integration of the dynamical system given by Eq. ^ 
in scale-free networks. 



IV. CONCLUSIONS 



pie models aiming at the reconstruction of the anatomical 
network based on the observed neurons' states will link, 
for this dynamics, neurons with similar in-degree with no 
connection in the real network (see Fig. [2]). 

The heterogeneous mean field solution allows us to re- 
late the degree distributions of the structural and func- 
tional topologies. In the remainder we will assume undi- 
rected structural networks. In the functional network a 
node with degree k is connected with a node with de- 
gree k' if their phase difference is smaller than a given 
threshold: |$fe - $fe'| < A. Then, using Eq. ([H]), the 
functional degree of a node of structural degree k is given 
by (l{k) ~ /|i/fe_i/fc'|<5 ^(^')^^' where 6 is an arbitrary 
threshold 5uj{k)T — A(l -I- {k)T)). If the degree distribu- 
tion of the structural network is a power law P{k) ^ k~^ 
then, q{k) ^ k^ where /3 = 2 — 7for7<2 and the degree 
distribution of the functional network is also power law 
P{q) ~ where a = (2 — 7)""'^. The numerical simula- 
tions of the system given by Eq. ([IJ shows an excellent 



Summarizing, we have got insight in the relationship 
between the topological network of connections and the 
functional network obtained from a simple dynamical 
process with delays. We have found the conditions for 
the emergence of locked dynamical states in any net- 
work of diffusively delay-coupled oscillators. We iden- 
tify these states as the main components of the emergent 
functional network generated by this simplified dynam- 
ics. Using these analytical guides we have explored the 
functional network obtained for the class of uncorrelated 
heterogenous networks, under the mean-field hypothesis, 
and have checked its prediction in scale-free networks. 
The results allow us to grasp the dependence of the func- 
tional network on the topological parameters, highlight- 
ing the role played by the delays and heterogeneity 35^]. 
Indeed, although functional and structural topologies dif- 
fer at the local level, we have shown that the degree dis- 
tributions are related in the presence of delays as distant 
nodes sharing the same degree will be functionally corre- 
lated. 



5 



V. ACKNOWLEDGMENTS 

V.M.E., A. A. and J.B.-H. acknowledge financial sup- 
port from MEC (Spain) through projects FISICOS 
(FIS2007-60327) and FIS2009-13730-C02-02. T.P. ac- 



knowledges support from the NSF, and Mathers Foun- 
dation. A. A. acknowledges partial support by the Direc- 
tor, Office of Science, Computational and Technology Re- 
search, U.S. Department of Energy under Contract No. 
DE-AC02-05CH11231. 



[1] A. Pikovsky, M. Rosenblum, and J. Kurtlis, Synchroniza- 
tion: A universal concept in nonlinear sciences (Cam- 
bridge University Press, Cambridge, England, 2003). 

[2] S. Strogatz, Sync: The emerging science of spontaneous 
order (Hyperion, New York, 2003). 

[3] A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno, and 
C. Zhou, Phys. Rep. 469, 93 (2008). 

[4] J.M. Stuart, E. Segal, D. Roller, SK Kim, Science, 302, 
249 (2003). 

[5] L. Mao, J.L. Van Hemert, S. Dash, and J. A. Dickerson, 

BMC Bioinformatics 10, 346 (2009). 
[6] A.G. Siapas, E.V. Lubenov, and M.A. Wilson, Neuron 

46, 141 (2005). 
[7] CM. Gonzalez, C. MasoUer, M.C. Torrent, and J. 

Garci'a-Ojalvo, Europhys. Lett. 79, 64003 (2007). 
[8] T. Gregor, K. Fujimoto, N. Masaki, and S. Sawai, Science 

328, 1021 (2010). 
[9] T. Danino, O. Mondragon-Palomino, L. Tsimirin, and J. 

Hasty, Nature 462, 326 (2010). 
[10] M. Chavez, M. Valencia, V. Navarro, V. Latora, and J. 

Martinerie, Phys. Rev. Lett. 104, 118701 (2010). 
[11] E. BuUmore and O. Sporns, Nat. Rev. Neurosci. 10, 186 

(2009) . 

[12] V.M. Egm'luz, D.R. Chialvo, G.A. Cecchi, M. Baliki, and 

A.V. Apkarian, Phys. Rev. Lett. 94, 018102 (2005). 
[13] M. Dhamala, V.K. Jirsa, and M.Z. Ding, Phys. Rev. Lett. 

92, 074104. (2004). 
[14] T. Y.-C. Tsai, Y. S. Choi, W. Ma, J. R. Pomerening, 

C. Tang, and James E. FerreU, Jr. Science 321, 126-129 

(2008). 

[15] L.G. Morelli, S. Ares, L. Herrgen, C.Schroter, F. Julicher, 
and A. C. Gates, The HFSP J. 3, 55-66 (2009). 

[16] L. Herrgen, S. Ares, L.G. MoreUi, C. Schroter, F. 
Jiilicher, and A.G. Gates, Current Biol. 20, 1244-1253 

(2010) . 

[17] D. Yu and S. Boccaletti, Phys. Rev. E 80, 036203 (2009). 
[18] I. Fischer, R. Vicente, J.M. Buldu, M. Peil, C.R. Mirasso, 

M.C. Torrent and J. Garcia-Ojalvo, Phys. Rev. Lett. 97, 

123902 (2006). 

[19] R. Vicente, L.L. GoUo, C.R. Mirasso, I. Fischer and P. 
Gordon, Proc. Natl. Acad. Sci. USA 105, 17157 (2008). 



[20] W. Kinzel, A. Englert, G. Reents, M. Zigzag and L Ran- 
ter, Phys Rev E 79, 056207 (2009). 

[21] G.L. Gerstein and A.M. Aertsen, J. Neurophysiol. 54 
1513 (1985). 

[22] G. Palm, A. Aertsen and G.L. Gerstein Biol. Cybern. 59 
1 (1988). 

[23] C. Zhou, L. Zemanova, G. Zamora, C.C. Hilgetag, J. 

Kurths, Phys. Rev. Lett. 97 238103 (2006). 
[24] D. Yu, M. Righero and L. Kocarev, Phys. Rev. Lett. 97, 

188701 (2006). 
[25] M. Timme, Phys. Rev. Lett. 98, 224101 (2007). 
[26] L Pajevic and D. Plenz PLoS Comp. Biol. 5, el000271 

(2008). 

[27] D. Napoletani and T. D. Sauer, Phys. Rev. E 77, 026103 

(2008) . 

[28] B. Barzel and O. Biham, Phys. Rev. E. 80, 046104 

(2009) . 

[29] Y. Kuramoto, in Lectures notes in Physics No 30, edited 

by H. Araki (Springer, New York, 1975). 
[30] F.R.K. Chung, CBMS Regional Conference Series in 

Mathematics 92, 1 (1997). 
[31] R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. 

Chklovskii and U. Alon, Science 298, 824 (2002). 
[32] O. Mason and M. Verwoerd, Systems Biology lET 1, 89 

(2007). 

[33] Y. Moreno, M. Vazquez-Prada, and A.F. Pacheco, Phys- 
ica A 343, 279 (2004). 

[34] Methods designed to extract structural connectivity have 
been proposed, for example, based on partial phase syn- 
chronization [B. Schelter, M. Winterhalder, R. Dhalhaus, 
J. Kurths, and J. Timmer, Phys. Rev. Lett. 96, 208103 
(2006)], on response dynamics to external driving 25], or 
on the presence of noise [J. Ren, W.-X. Wang, B. Li, and 
Y.-C. Lai, Phys. Rev. Lett. 104, 058701 (2010)]. 

[35] The relation obtained between the functional and struc- 
tural degree distribution suggests an interesting link: 
functional brain networks have been observed to display 
exponents close to 2 [l2], with anatomical neuronal net- 
works showing an exponent close to 1.6 [G. Gong et al, 
Cereb. Cortex 19, 524-536 (2009)]. 



